setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file2=read.csv("./20210120.H3k.analysis/117012.analysis/117K.ASH.add.enh.promtr.csv",header=T)
names(file2)=c("unitID","snp.location")
filea=read.csv("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)

asm2=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/220520ASMs_anno.hg19_multianno.csv",head=T,sep=",")	#220K
asm2$unitID=paste(asm2$Chr,asm2$Start,asm2$Ref,asm2$Alt,sep=":")
asm=read.csv("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/869727.all.snp.vaf.up.down.20210321.csv",head=T)
asm=asm[asm$unitID %in% asm2$unitID,]
asm=asm[asm$unitID %in% filea$unitID,]#13649行
asm=merge(file2,asm,by="unitID")
asm=asm[abs(asm$chazhi)>0.1,]
asm=data.frame(unitID=asm$unitID,ASM.alt.group=asm$ASM.alt.group,snp.location=asm$snp.location)

file1=read.table("./20210321.H3K分析/117012.vaf.up.down.20210321.txt",header=T,sep="\t")
file1=merge(file1,file2,by="unitID")
file1=file1[abs(file1$chazhi)>0.1,]
file1=data.frame(unitID=file1$unitID,ASH.alt.group=file1$alt.group,snp.location=file1$snp.location)

rt=merge(file1[,1:2],asm,by="unitID")#去掉其中一项snp.location
rt$direction=paste(rt$ASH.alt.group,rt$ASM.alt.group,sep=":")
group1=c("enhancer","promoter","")

for(i in 1:3){
tmp=rt[rt$snp.location==group1[i],]
result=data.frame(matrix(,nrow = 2,ncol = 2))
names(result)=c("ASM.up","ASM.down")
row.names(result)=c("ASH.down","ASH.up")
result[1,1]=dim(tmp[tmp$direction=="down:up",])[1]
result[1,2]=dim(tmp[tmp$direction=="down:down",])[1]
result[2,1]=dim(tmp[tmp$direction=="up:up",])[1]
result[2,2]=dim(tmp[tmp$direction=="up:down",])[1]

fn=paste("./20210321.H3K分析/位点羟甲基化和甲基化对比.",group1[i],".csv",sep="")

write.csv(result,fn,,quote=F,row.names=T)
}
